## ---- meta
if (.Platform$OS.type == "windows") {
  Sys.setenv(ADMB_HOME="C:\\progra~1\\ADMB",
             BCC55_HOME="C:\\Borland\\BCC55",
             PATH=paste("C:\\Borland\\BCC55\\bin;C:\\progra~1\\ADMB\\bin;", shell("echo %PATH%", intern = TRUE), sep = ""))
}

source("../../code/R/functions_OpModel.R")
source("../../code/R/functions_simulate.R")
source("../../code/R/functions_ADMB.R")

system("simulate_RE.tpl")
system("estimate_nlme.tpl")
system("estimate_bay.tpl")

#system("cp ../../code/ADMB/simulate_RE.tpl .")
#system("cp ../../code/ADMB/estimate_nlme.tpl .")
#system("cp ../../code/ADMB/estimate_bay.tpl .")

save(list = ls(all = TRUE), file = "100simsMuPlus1DepPt4sigRpt1_lag.Rda")

## ---- parameter determination
dep <- .2 # depletion relative to B0
reps <- 20
T_dep <- 20

mu <- 0
tau <- 1

T <- 40
m <- 10
beta <- rnorm(m, mu, tau)
h <- (exp(beta) + 0.2) / (1 + exp(beta))

Eb <- rep(c(rep(0, len=T/8), seq(0,.1, len=T/8), seq(.1,.25, len=T/4), seq(.25, 0, len=T/4), rep(0, len=T/4), 0))
E <- matrix(rep(Eb, times = m), T+1, m, byrow = F)

L <- 4 # try making this diff for each stock as well
B0 <- rep(40, m); p <- rep(1, m); w <- rep(1, m); M <- rep(0.2, m);
sigR <- 0.4; qI <- 1; sigI <- 0.5; q <- rep(1, m); sigN <- 0.5; qN <- 1

phi0 <- rep(0,m) # mle estimates
for (k in 1:m)
  phi0[k] <- (1 - p[k]*w[k]*exp(-M[k])) /
  (1 - (1+p[k])*exp(-M[k]) + p[k] * exp(-2*M[k]))

R0 <- B0 / phi0

pars <- list(beta=beta, mu=mu, tau=tau, B0=B0, p=p, w=w, M=M, sigR=sigR,
             qI=qI, sigI=sigI, sigN=sigN, qN=qN, q=q)

sim <- opmodel(pars, m, T, E, L, phi0)
errR <- sim$errR; errI <- sim$errI; errN <- sim$errN

## ---- loop around and optimize q, simulate!
q.temp <- rep(0, length = m)
for (k in 1:m) {
  q.temp[k] <- 
    optimize(opmodel.dep, interval = c(.1,50), T_dep = T_dep, dep = dep, k = k,
             pars = pars, m = m, T = T, E = E, L = L, phi0 = phi0,
             errR = errR, errI = errI, errN = errN)$minimum
}

pars$q <- q.temp
sim <- opmodel(pars, m, T, E, L, phi0, errR, errI, errN)

meanR0 <- apply(sim$R * apply(sim$B, 2, function(x) x > median(x)), 2, function(x) mean(x[x>0]))

## do the meta-analysis methods (NLME)
system("admb -r simulate_RE"); system("admb -r estimate_nlme");

results <- list()
sims <- 10
for (i in 1:sims) {
  pars$beta <- rnorm(m, mu, tau)
  sim <- run_opmodel(pars, m, T, E, L, phi0)
  #est <- run_estmodel(sim, pars)
  out <- run_meta(est, pars, method = "nlme", real = TRUE)
  results[[i]] <- data.frame(beta = pars$beta,
                             beta_est = with(out, eps * exp(log_tau) + mu),
                             eps = out$eps,
                             B0 = B0,
                             B0_est = out$B0,
                             mu = mu,
                             mu_est = out$mu,
                             tau = tau,
                             tau_est= exp(out$log_tau),
                             sigR = sigR,
                             #sigR_sim = as.numeric(lapply(est, function(x) x$sigR)),
                             sigR_est = exp(out$log_sigR))
  cat(paste("\n", i, "asdkfj;askldjfa;lksdf;askdf;\n"))
}

## plot performances of simulations of all stocks
aho <- do.call(rbind, results)

par(mfrow = c(3,3), mar = c(5-1,4,4-1,2))
plot(beta_to_h(aho$beta), beta_to_h(aho$beta_est), xlim = c(.2,1), ylim = c(.2,1), 
     col = rgb(0,0,0,.3), pch = 16, main = "steepness", cex = 2, xlab = "actual", ylab = "predicted")
abline(a = 0, b = 1)
hist(aho$B0_est, breaks = 30, main = "B0", col = "grey", border = "lightgrey")
abline(v = unique(aho$B0))
hist(aho$sigR_est, breaks = 30, main = "sigR", col = "grey", border = "lightgrey")
abline(v = unique(aho$sigR))
dotchart(unique(aho$mu_est), main = "mu", xlim = c(-3,3)); abline(v = mu)
dotchart(unique(aho$tau_est), main = "tau", xlim = c(0,2)); abline(v = tau)
boxplot(aho$eps ~ as.character(rep(1:m, each = sims)), horizontal = TRUE, main = "eps")
plot(beta_to_h(aho$beta_est), aho$B0_est, col = rgb(0,0,0,.3), pch = 16, main = "", cex = 2)
plot(beta_to_h(aho$beta_est), aho$sigR_est, col = rgb(0,0,0,.3), pch = 16, main = "", cex = 2)
matplot(sim$B, type = "l", ylim = c(0,max(sim$B)))

## do the meta-analysis methods (NLME)
system("admb estimate_bay");

results <- list()
sims <- 10
for (i in 1:sims) {
  pars$beta <- rnorm(m, mu, tau)
  sim <- run_opmodel(pars, m, T, E, L, phi0)
  meanR0 <- apply(sim$R * apply(sim$B, 2, function(x) x > median(x)), 2, function(x) mean(x[x>0]))
  #est <- run_estmodel(sim, pars)
  out <- run_meta(est, pars, method = "bay", real = TRUE)
  results[[i]] <- data.frame(beta = pars$beta,
                             beta_mc = apply(out$beta_mc, 2, median),
                             B0 = B0,
                             B0_mc = apply(out$B0_mc, 2, median),
                             mu = mu,
                             mu_mc = median(out$mu_mc),
                             tau = tau,
                             tau_mc = median(exp(out$log_tau_mc)),
                             sigR = sigR,
                             #sigR_sim = as.numeric(lapply(est, function(x) x$sigR)),
                             sigR_mc = apply(exp(out$log_sigR_mc), 2, median))
  cat(paste("\n", i, "asdkfj;askldjfa;lksdf;askdf;\n"))
}

# for bay
aho <- do.call(rbind, results)

par(mfrow = c(3,3))
plot(beta_to_h(aho$beta), beta_to_h(aho$beta_mc), xlim = c(.2,1), ylim = c(.2,1), 
     col = rgb(0,0,0,.3), pch = 16, main = "steepness", cex = 2, xlab = "actual", ylab = "predicted")
abline(a = 0, b = 1)
hist(aho$B0_mc, breaks = 30, main = "B0", col = "grey", border = "lightgrey")
abline(v = unique(aho$B0))
hist(aho$sigR_mc, breaks = 30, main = "sigR", col = "grey", border = "lightgrey")
abline(v = unique(aho$sigR))
dotchart(unique(aho$mu_mc), main = "mu", xlim = c(-3,3)); abline(v = mu)
dotchart(unique(aho$tau_mc), main = "tau", xlim = c(0,2)); abline(v = tau)
plot(beta_to_h(aho$beta_mc), aho$B0_mc, col = rgb(0,0,0,.3), pch = 16, main = "", cex = 2)
plot(beta_to_h(aho$beta_mc), aho$sigR_mc, col = rgb(0,0,0,.3), pch = 16, main = "", cex = 2)
matplot(sim$B, type = "l", ylim = c(0,max(sim$B)))
